home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dgghrd.f < prev    next >
Text File  |  1996-11-04  |  8KB  |  254 lines
  1.       SUBROUTINE DGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
  2.      $                   LDQ, Z, LDZ, INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          COMPQ, COMPZ
  11.       INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
  15.      $                   Z( LDZ, * )
  16. *     ..
  17. *
  18. *  Purpose
  19. *  =======
  20. *
  21. *  DGGHRD reduces a pair of real matrices (A,B) to generalized upper
  22. *  Hessenberg form using orthogonal transformations, where A is a
  23. *  general matrix and B is upper triangular:  Q' * A * Z = H and
  24. *  Q' * B * Z = T, where H is upper Hessenberg, T is upper triangular,
  25. *  and Q and Z are orthogonal, and ' means transpose.
  26. *
  27. *  The orthogonal matrices Q and Z are determined as products of Givens
  28. *  rotations.  They may either be formed explicitly, or they may be
  29. *  postmultiplied into input matrices Q1 and Z1, so that
  30. *
  31. *       Q1 * A * Z1' = (Q1*Q) * H * (Z1*Z)'
  32. *       Q1 * B * Z1' = (Q1*Q) * T * (Z1*Z)'
  33. *
  34. *  Arguments
  35. *  =========
  36. *
  37. *  COMPQ   (input) CHARACTER*1
  38. *          = 'N': do not compute Q;
  39. *          = 'I': Q is initialized to the unit matrix, and the
  40. *                 orthogonal matrix Q is returned;
  41. *          = 'V': Q must contain an orthogonal matrix Q1 on entry,
  42. *                 and the product Q1*Q is returned.
  43. *
  44. *  COMPZ   (input) CHARACTER*1
  45. *          = 'N': do not compute Z;
  46. *          = 'I': Z is initialized to the unit matrix, and the
  47. *                 orthogonal matrix Z is returned;
  48. *          = 'V': Z must contain an orthogonal matrix Z1 on entry,
  49. *                 and the product Z1*Z is returned.
  50. *
  51. *  N       (input) INTEGER
  52. *          The order of the matrices A and B.  N >= 0.
  53. *
  54. *  ILO     (input) INTEGER
  55. *  IHI     (input) INTEGER
  56. *          It is assumed that A is already upper triangular in rows and
  57. *          columns 1:ILO-1 and IHI+1:N.  ILO and IHI are normally set
  58. *          by a previous call to DGGBAL; otherwise they should be set
  59. *          to 1 and N respectively.
  60. *          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
  61. *
  62. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
  63. *          On entry, the N-by-N general matrix to be reduced.
  64. *          On exit, the upper triangle and the first subdiagonal of A
  65. *          are overwritten with the upper Hessenberg matrix H, and the
  66. *          rest is set to zero.
  67. *
  68. *  LDA     (input) INTEGER
  69. *          The leading dimension of the array A.  LDA >= max(1,N).
  70. *
  71. *  B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)
  72. *          On entry, the N-by-N upper triangular matrix B.
  73. *          On exit, the upper triangular matrix T = Q' B Z.  The
  74. *          elements below the diagonal are set to zero.
  75. *
  76. *  LDB     (input) INTEGER
  77. *          The leading dimension of the array B.  LDB >= max(1,N).
  78. *
  79. *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
  80. *          If COMPQ='N':  Q is not referenced.
  81. *          If COMPQ='I':  on entry, Q need not be set, and on exit it
  82. *                         contains the orthogonal matrix Q, where Q'
  83. *                         is the product of the Givens transformations
  84. *                         which are applied to A and B on the left.
  85. *          If COMPQ='V':  on entry, Q must contain an orthogonal matrix
  86. *                         Q1, and on exit this is overwritten by Q1*Q.
  87. *
  88. *  LDQ     (input) INTEGER
  89. *          The leading dimension of the array Q.
  90. *          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
  91. *
  92. *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
  93. *          If COMPZ='N':  Z is not referenced.
  94. *          If COMPZ='I':  on entry, Z need not be set, and on exit it
  95. *                         contains the orthogonal matrix Z, which is
  96. *                         the product of the Givens transformations
  97. *                         which are applied to A and B on the right.
  98. *          If COMPZ='V':  on entry, Z must contain an orthogonal matrix
  99. *                         Z1, and on exit this is overwritten by Z1*Z.
  100. *
  101. *  LDZ     (input) INTEGER
  102. *          The leading dimension of the array Z.
  103. *          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
  104. *
  105. *  INFO    (output) INTEGER
  106. *          = 0:  successful exit.
  107. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  108. *
  109. *  Further Details
  110. *  ===============
  111. *
  112. *  This routine reduces A to Hessenberg and B to triangular form by
  113. *  an unblocked reduction, as described in _Matrix_Computations_,
  114. *  by Golub and Van Loan (Johns Hopkins Press.)
  115. *
  116. *  =====================================================================
  117. *
  118. *     .. Parameters ..
  119.       DOUBLE PRECISION   ONE, ZERO
  120.       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  121. *     ..
  122. *     .. Local Scalars ..
  123.       LOGICAL            ILQ, ILZ
  124.       INTEGER            ICOMPQ, ICOMPZ, JCOL, JROW
  125.       DOUBLE PRECISION   C, S, TEMP
  126. *     ..
  127. *     .. External Functions ..
  128.       LOGICAL            LSAME
  129.       EXTERNAL           LSAME
  130. *     ..
  131. *     .. External Subroutines ..
  132.       EXTERNAL           DLARTG, DLASET, DROT, XERBLA
  133. *     ..
  134. *     .. Intrinsic Functions ..
  135.       INTRINSIC          MAX
  136. *     ..
  137. *     .. Executable Statements ..
  138. *
  139. *     Decode COMPQ
  140. *
  141.       IF( LSAME( COMPQ, 'N' ) ) THEN
  142.          ILQ = .FALSE.
  143.          ICOMPQ = 1
  144.       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
  145.          ILQ = .TRUE.
  146.          ICOMPQ = 2
  147.       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
  148.          ILQ = .TRUE.
  149.          ICOMPQ = 3
  150.       ELSE
  151.          ICOMPQ = 0
  152.       END IF
  153. *
  154. *     Decode COMPZ
  155. *
  156.       IF( LSAME( COMPZ, 'N' ) ) THEN
  157.          ILZ = .FALSE.
  158.          ICOMPZ = 1
  159.       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
  160.          ILZ = .TRUE.
  161.          ICOMPZ = 2
  162.       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
  163.          ILZ = .TRUE.
  164.          ICOMPZ = 3
  165.       ELSE
  166.          ICOMPZ = 0
  167.       END IF
  168. *
  169. *     Test the input parameters.
  170. *
  171.       INFO = 0
  172.       IF( ICOMPQ.LE.0 ) THEN
  173.          INFO = -1
  174.       ELSE IF( ICOMPZ.LE.0 ) THEN
  175.          INFO = -2
  176.       ELSE IF( N.LT.0 ) THEN
  177.          INFO = -3
  178.       ELSE IF( ILO.LT.1 ) THEN
  179.          INFO = -4
  180.       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
  181.          INFO = -5
  182.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  183.          INFO = -7
  184.       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  185.          INFO = -9
  186.       ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
  187.          INFO = -11
  188.       ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
  189.          INFO = -13
  190.       END IF
  191.       IF( INFO.NE.0 ) THEN
  192.          CALL XERBLA( 'DGGHRD', -INFO )
  193.          RETURN
  194.       END IF
  195. *
  196. *     Initialize Q and Z if desired.
  197. *
  198.       IF( ICOMPQ.EQ.3 )
  199.      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
  200.       IF( ICOMPZ.EQ.3 )
  201.      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
  202. *
  203. *     Quick return if possible
  204. *
  205.       IF( N.LE.1 )
  206.      $   RETURN
  207. *
  208. *     Zero out lower triangle of B
  209. *
  210.       DO 20 JCOL = 1, N - 1
  211.          DO 10 JROW = JCOL + 1, N
  212.             B( JROW, JCOL ) = ZERO
  213.    10    CONTINUE
  214.    20 CONTINUE
  215. *
  216. *     Reduce A and B
  217. *
  218.       DO 40 JCOL = ILO, IHI - 2
  219. *
  220.          DO 30 JROW = IHI, JCOL + 2, -1
  221. *
  222. *           Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
  223. *
  224.             TEMP = A( JROW-1, JCOL )
  225.             CALL DLARTG( TEMP, A( JROW, JCOL ), C, S,
  226.      $                   A( JROW-1, JCOL ) )
  227.             A( JROW, JCOL ) = ZERO
  228.             CALL DROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
  229.      $                 A( JROW, JCOL+1 ), LDA, C, S )
  230.             CALL DROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
  231.      $                 B( JROW, JROW-1 ), LDB, C, S )
  232.             IF( ILQ )
  233.      $         CALL DROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C, S )
  234. *
  235. *           Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
  236. *
  237.             TEMP = B( JROW, JROW )
  238.             CALL DLARTG( TEMP, B( JROW, JROW-1 ), C, S,
  239.      $                   B( JROW, JROW ) )
  240.             B( JROW, JROW-1 ) = ZERO
  241.             CALL DROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
  242.             CALL DROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
  243.      $                 S )
  244.             IF( ILZ )
  245.      $         CALL DROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )
  246.    30    CONTINUE
  247.    40 CONTINUE
  248. *
  249.       RETURN
  250. *
  251. *     End of DGGHRD
  252. *
  253.       END
  254.